Diffusive Transport Enhanced by Thermal Velocity Fluctuations 
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We study the contribution of advection by thermal velocity fluctuations to the effective diffusion 
coefficient in a mixture of two indistinguishable fluids. We find good agreement between a simple 
fluctuating hydrodynamics theory and particle and finite-volume simulations. The enhancement of 
the diffusive transport depends on the system size L and grows as ln(L/Lo) in quasi two-dimensional 
systems, while in three dimensions it scales as Lq 1 — L^ 1 , where Lq is a reference length. Our results 
demonstrate that fluctuations play an important role in the hydrodynamics of small-scale systems. 



Thermal fluctuations in non-equilibrium systems in 
which a constant (temperature, concentration, velocity) 
gradient is imposed externally exhibit remarkable be- 
havior compared to equilibrium systems jlj. The so- 
lution of the linearized equations of fluctuating hydro- 
dynamics shows that concentration and density fluctu- 
ations exhibit long-ranged correlations in the presence 
of a macroscopic concentration gradient Vc 1-3 . The 
enhancement of large-scale (small wavenumber) concen- 
tration fluctuations is dramatic during the early stages of 
diffusive mixing between initially phase-separated fluids. 
These giant fluctuations ^[5] during free diffusive mix- 
ing have been measured using light scattering and shad- 
owgraphy techniques [H [7] , finding good but imperfect 
agreement with theoretical predictions. 

The giant fluctuation phenomenon arises because of 
the appearance of long-ranged correlations between con- 
centration and velocity fluctuations in the presence of a 
concentration gradient. It has been predicted that these 
correlations give rise to fluctuation-renormalized trans- 
port coefficients [3l[8]; however, the predicted enhance- 
ment of transport at hydrodynamic scales has not yet 
been computationally observed. In particular, it is im- 
portant to understand how the effective transport coeffi- 
cients depend on the length scale of observation. 

In this Letter we consider diffusion in a mixture of 
identical but labeled (as components 1 and 2) fluids 
[9] enclosed in a box of size L x x L y x L z , in the ab- 
sence of gravity. Periodic boundary conditions are ap- 
plied in the x (horizontal) and z (depth) directions, 
while the top and bottom boundaries are impermeable 
constant-temperature walls. A concentration gradient 
Vc = {cT — CB)/L y is imposed along the y axes by enforc- 
ing a constant concentration ct at the top wall and Cb 
at the bottom wall. Because the fluids are indistinguish- 
able, concentration is passively transported by thermal 
fluctuations. 

Since species are not changed in particle collisions, the 
diffusive transport of concentration c = pi/p can only 
occur via advective motion of the particles, where p de- 
notes the mass density. The mass flux for a given species 
is therefore equal to the momentum density for particles 



of that species. At steady state the particles of a given 
species have a non-zero macroscopic momentum density 
ji = p\Vx = — px(Vc), where x is t ne mass diffusion 
coefficient [TIT] • The local fluctuations around the macro- 
scopic mean, p\ = p\ + 6p\ and V\ — V\ + 8v\, can also 
make a non-trivial contribution to the average mass flux 
if they are correlated, 

<Ji) = {pxvi) = -PX (Vc) + ((S Pl ) (S Vl )) . (1) 
At mesoscopic scales the hydrodynamic behavior of flu- 
ids can be described with the Landau-Lifshitz Navier- 
Stokes (LLNS) equations of fluctuating hydrodynamics 
[TJ [TT] . The incompressible isothermal LLNS equations 
for a mixture of two indistinguishable fluids are 

d t v= -Vir -v-Vv + i>V 2 v + V-(A V W), (2) 

d t c = -v-Vc + x 'V 2 c + V-(A C \V) (3) 
where rj is the viscosity and v = rj/p, and the pressure it 
enforces 'V-v = 0. The stochastic fluxes are white-noise 
random Gaussian tensor W and vector W fields, with 
amplitudes A 2 = 2nksT '/ p 2 and A 2 = 2mxc(l — c)/p 
determined from the fluctuation-dissipation principle, 
where m is the fluid particle mass. 

In addition to the usual Fickian contribution, the dif- 
fusive flux in (|3| includes advection by the fluctuating 
velocities v = Sv, 

-v • V (6c) + xV 2 (6c) = V ■ [- (6c) (Sv) + xV (6c)} , 
which is a quadratic function of the fluctuations. To lead- 
ing order, the advective contribution to the average dif- 
fusive mass flux is approximated using the solution of the 
linearized equations, 

-((6c) (Sv)) « -((6c) (6v)) linca r = (A X ) Vc. 
The effective diffusion coefficient XcK = X + Ax thus in- 
cludes an enhancement Ax due to thermal velocity fluc- 
tuations, in addition to the hare diffusion coefficient x- 

The solution to the linearized form of ( 2j3 ) in the 
Fourier domain [3J [TJ] shows that the concentration fluc- 
tuations and the fluctuations of velocity parallel to the 
gradient develop long ranged correlations, 

S c , n = ((6c)(Tvl)) = - J B + T x)k2 (sm 2 9) Vc. (4) 

where 9 is the angle between k and Vc, sin 2 = k\jk 2 , 
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a hat denotes the Fourier transform, and star denotes 
the complex conjugate. The power-law divergence for 
small k indicates long ranged correlations between 8c and 
5vu, and is the cause of both the giant fluctuation phe- 
nomenon and the diffusion enhancement. As seen from 
pi), the actual correlation that determines the diffusion 

. ~ (!) 

enhancement is S p — ((5c)(Sv^ )*) s» pS CiV ^. 

We verify the predictions of fluctuating hydrodynam- 
ics by using the Direct Simulation Monte Carlo (DSMC) 
particle algorithm |13j . Previous careful measurements 
of transport coefficients in DSMC have been limited to 
quasi one-dimensional simulations |14j . The effect we are 
exploring here does not appear in one dimension as it 
arises because of the presence of vortical modes in the 
fluctuating velocities. We have performed DSMC calcu- 
lations for an ideal hard-sphere gas with molecular di- 
ameter cr = l and molecular mass m = 1, at an equi- 
librium density of p = 0.06, with the temperature kept 
at fcsTo = To = 1 via thermal collisions with the top 
and bottom walls. Each DSMC particle represents a sin- 
gle hard sphere so the mean free path is A = 3.75 and 
the mean free collision time is r = 2.35. The DSMC time 
step was chosen to be At = t/2, and the collision cell size 
is either Ax c — A or Ax c = 2A. A uniform concentration 
gradient along the vertical (y) direction is implemented 
by randomly selecting the species of particles to be one 
with probability ct/b if they collide with the top/bottom 
wall. Hydrodynamic quantities such as velocity and con- 
centration are calculated from the particle data by using 
a grid of N x x N y x N z sampling or hydrodynamic cells, 
each of volume AV = Ax Ay Az, and a discrete Fourier 
transform is used to obtain static structure factors. 

To compare the prediction Q to results from parti- 
cle simulations, we have converted the continuum static 
structure factor iS c ,U|i (&) hito a discrete structure fac- 
tor Sc,v^ for finite-volume averages of the continuum 
fields, where the wavenumbers k£Z 3 index the discrete 
set of wavevectors compatible with periodicity [12] . In 
Fig. fljwe compare the theoretical prediction for pS CyV ^ (k) 
to DSMC results for the discrete structure factor S m . 

Similar results are obtained for two different sizes of the 
DSMC collision cells [12], Ax c = 2A and Ax c = A, verify- 
ing that the details of the microscopic collision dynamics 
do not affect the mesoscopic hydrodynamic behavior. 

It is expected that compressibility effects would affect 
S (i). In order to construct a theoretical prediction, 

however, one must not only include the effect of com- 
pressibility but also replace the "one-fluid" approxima- 
tion with a corresponding "two-fluid" compressible hy- 
drodynamic theory [10]. As Fig. [I] demonstrates, the 
incompressible isothermal theory for pS CjV ,, can be used 
as a proxy for S m in order to construct theoretical 

predictions for the diffusion enhancement. 

The mass flux due to advection by the fluctuating ve- 
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Figure 1: (Color online) Discrete structure factor 5 <i) from 

quasi two-dimensional DSMC runs with L x — 64A, L y = 512A 
and L z = 2A, for several wavenumbers km = K x ■ 2n / L x (see 
legend) , compared to the discrete equivalent of the continuum 
prediction Q (solid lines of the same color). Note that for a 
fixed k x we expect the structure factor to decay as ky 4 . 



locities can be approximated using Q as 
((5c) (Sv)) 

linear — 

(27r) _J / S C) „ (k) dk, 

Jk 



(5) 



giving an estimate of the diffusion enhancement [3) I12j 

A X = , , 3 fc f\ , / ( si * 2 0) k ^ dk- (6) 
(2tt)'V (X + v) J k 

Because of the fc _2 -like behavior, the integral over all 

k above diverges unless one imposes [3] a lower bound, 

fcmin ~ 2ir/L in the absence of gravity, and a phenomeno- 

logical cutoff fc max ~ 7r/L mo i for the upper bound, where 

L mo i is an ad-hoc "molecular" length scale. 

For a quasi two-dimensional system, L z < Lj «L 9 , 
we can replace the integral over k z with 2ir/L z and inte- 
grate over all k y . This leads to an average total diffusive 
flux that grows logarithmically with the width L x for a 
fixed height L y [12], 



(2D) 
Xeff 



X+ 4Mx + v)L z ^ (7) 
where Lq > 2L mo i is the reference width at which 
the "bare" diffusion coefficient is measured. Within 
the phenomenological perturbative theory Lq is an ar- 
bitrary (mesoscopic) length scale, and simply defines 
X = Xcff(L x = ^o)- For comparison between the particle 
simulations and the theory we use Lq = 16A. When the 
system width becomes comparable to the height, bound- 
aries will intervene and for L x 3> L y the effective dif- 
fusion coefficient must become a constant, which is pre- 
dicted to be a logarithmically-growing function of L y in 
two dimensions. The same coefficient in front of lni^ as 
in ([7| is obtained when the integral over k x is replaced 
by a discrete sum over the wavenumbers consistent with 
periodicity, k x = k x ■ 2tt/L x , k x 



In three dimensions, L 7 



L z — L < L y , Xeff con- 
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(8) 



verges as L — > oo to the macroscopic diffusion coefficient, 
(3D) ak B T / 1 1 

Xcff WX %(x + ^)Uo 
but for a finite system the effective diffusion coefficient 
is reduced by an amount ~ L^ 1 due to the truncation 
of the velocity fluctuations by the confining walls. Cal- 
culating the exact value of a requires performing a sum 
over k x and k z instead of integrals over k x and k z , as we 
have done numerically |f 2j . The numerical results sug- 
gest that, as in two dimensions, the difference in xiff°^ 
between two systems attains a finite value as L mo \ — > 0, 
justifying Q for (L ,L) > L mol . 

In particle simulations, we calculate the effective dif- 
fusion coefficient Xeff from the momentum density of one 
of the species (denoted either with a subscript or with a 
parenthesis superscript) along the vertical direction, 



PoXcS 



cr - c B 



Li, 



Ay 



PoXeff(Vc), (9) 



where we measure ct and cb in the top and bottom layer 
of sampling cells (whose centers are a distance L y — Ay 
from each other) to empirically account for the small con- 
centration slip in DSMC. Numerical experiments have 
verified that ) matches the flux obtained from count- 
ing the average number of color flips at the top or bottom 
walls. Furthermore, the results verify that (in ) is linear 
in the gradient Vc, and that c~t/b ar e essentially inde- 
pendent of the system dimensions. 

The traditional definition of a "renormalized" diffusion 
coefficient [8] as the macroscopic limit of Xeffj only works 
in three dimensions and is not very useful for confined 
systems. Instead, for each sampling cell, we define a lo- 
cally renormalized diffusion coefficient xo via 

(pi)<«5 1) ) = <Pi)0{ 1) M>=PXo(Vc) ! (10) 

where we have accounted for the fact that the macro- 
scopic concentration profile c(y) may depend on y. In 
fact, such a dependence is observed in the particle simula- 
tions, and we have approximated the local concentration 
gradient dc/dy by a numerical derivative of a polynomial 
of degree five fit to c(y). We have empirically observed 
that xo is independent of y, except for a boundary layer 
close to the top and bottom walls [12]. This is an im- 



portant finding, since ( 10 ) is a constitutive model that 



is assumed to hold not just at the macroscale but also 
at the mesoscale, notably, xo is an input parameter for 
fluctuating hydrodynamics finite- volume solvers [15] . 

Figure [2^ shows how Xeff and xo change as the width 
of the system L x is increased while keeping the height L y 
fixed for two different quasi two-dimensional DSMC sys- 
tems. For System A, the DSMC collision cells are cubes 
of side Ax c = 7.5 = 2A, while each sampling cell contains 
2x2x1 collision cells, or N p — 101 particles on aver- 
age. The height of the box is L y — 256A = 960 and the 
imposed concentrations at the walls are cb = 0.25 and 
ct = 0.75. For System B, the sampling cells are twice 



as large, 4x4x1 collision cells each, and the system 
height is L y = 512A = 1920. We obtain similar results 
using a factor of two smaller collision cells (not shown). 
For the quasi two-dimensional systems, the thickness is 
L z = 7.5 = 2A and there is only one DSMC collision cell 
along the z direction. Figure[2^L shows that Xeff grows like 
\nL x , with a slope that is well-predicted by Eq. Q. For 
widths larger than about 8 mean free paths, xo becomes 
constant and rather similar to the Chapman-Enskog ki- 
netic theory prediction. Note that xo is not a fundamen- 
tal material constant and in fact depends on the shape 
of the sampling cells, notably, it grows as the sampling 
cell size is enlarged. 

In Fig. [2]d we show results from three dimensional 
DSMC simulations, in which the system width (%) and 
depth (z) directions are equivalent, L z = L x = L, and 
the rest of the parameters are as for System A. Similar 
behavior is seen as in two dimensions, except that now 
the effective diffusion grows as — L^ 1 and saturates to a 
constant value for large L, assuming that L y L. 

The predictions of the simplified fluctuating hydrody- 
namic theory, Eqs. ^ and (J8|, are shown in Figs. [2] and 
are seen to be in very good agreement with the particle 
simulations for intermediate L x . However, recall that the 
incompressible isothermal theory assumed that L y is es- 
sentially infinite and thus in two dimensions Xeff grows 
unbounded in the macroscopic limit. Yet when L x S> L y , 
Xeff must saturate to a constant value, and the parti- 
cle data shown in Fig. [2^, shows measurable deviations 
from the simple theory for L x > L y /2. One can extend 
the theoretical calculations to account for the hard wall 
boundary conditions in the y direction p], however, such 
a calculation is non trivial. Instead, we have used the 
finite-volume solver developed in Ref. [TS] to solve the 
LLNS non-linear system of SPDEs for the same system 
dimensions as in the particle simulations. To minimize 
the effect of nonlinearities in the SPDE solver, we arti- 
ficially reduce the amplitude of the noise by some factor 
e <C 1, but then scale all correlations by e~ 2 [12 . The 
results, shown in Fig. [2] are in excellent agreement with 
the particle simulations for the larger system sizes. 

In finite-volume solvers, the spacing of the computa- 
tional grid plays the equivalent of the cutoff length L mo i, 
and therefore Xeff depends on the grid spacing. We have 
added a constant to the effective diffusion coefficient ob- 
tained from SPDE runs so as to match Xeff from the par- 
ticle simulations for L x = Lq = 16A. This correction 
essentially renormalizes xo based on the size of the finite- 
volume hydrodynamic cells. One can think of xo as the 
physical-space equivalent of the wavenumber-dependent 
diffusion coefficient x(k,L<J = 0) commonly used in lin- 
ear response theories [H [9] . Theoretical predictions [12] 
for xo indicate that xo only includes "sub-grid" contribu- 
tions, from wavenumbers larger than 2ir/Ax. Thus xo 
stops increasing once the system becomes substantially 
larger than the size of the sampling cell. The bare diffu- 
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Figure 2: (Color online) (Panel a) The effective Xcff an d 
the renormalized xo diffusion coefficients as a function of the 
width of the system L x in two dimensions. Numerical results 
for System A (DSMC and SPDE) and System B (DSMC) are 
shown with symbols (see legend). The error bars for all of 
the simulation data are comparable or smaller than the size 
of the symbols. The theoretical predictions [12] are evaluated 
numerically and shown with lines. (Panel b) Same as panel 
(a) but in three dimensions. The inset highlights the L _1 
behavior. 

sion coefficient \ in the theory and SPDE calculations is 
adjusted so that for L x = Lq = 16A the effective diffusion 
is the same as that measured in the particle simulations. 

Previously-studied corrections to the bare or molecu- 
lar transport coefficients due to the tail of the velocity 
autocorrelation function (VACF) [2J, hydrodynamic in- 
teractions with periodic images of a given particle |16) . 
and the contribution due to advection by thermal ve- 
locity fluctuations 3, 8 J studied here, are all the same 
physical phenomenon simply calculated through differ- 
ent theoretical approaches, all of which are equivalent 
because of linearity |12j . In three dimensions, simple 
estimates indicate that the contribution of fluctuations 
to the macroscopic diffusion coefficient are small com- 
pared to molecular effects for gases but can be signif- 
icant for liquids [12]. However, the logarithmic diver- 
gence in (JtJ) means that A\ ^> X f° r sufficiently large 
(quasi) two-dimensional systems, requiring the inclusion 
of higher order corrections in the theory. At present, 



reaching the system width L x where A% ^ x is diffi- 
cult with DSMC simulations, but it may be accessible to 
finite- volume SPDE solvers or experiments [6]. 

Our results conclusively demonstrate that the advec- 
tion by thermal velocity fluctuations affects the mean 
transport in nonequilibrium finite systems. Theoretical 
modeling of finite systems at the nano or microscale thus 
requires including nonlinear hydrodynamic fluctuations. 
The advantage of fluctuating hydrodynamics is that it is 
simple, and it can take into account the proper boundary 
conditions and exact geometry, especially if a numerical 
SPDE solver is used. Furthermore, other effects such 
as gravity, temperature variations, or time dependence, 
can easily be included. However, a proper fully-nonlinear 
theory has yet to be developed, and requires detailed un- 
derstanding of the role of the necessary large wavenumber 
cutoffs (regularizations). Future work should verify the 
predictions of fluctuating hydrodynamics for the effect 
of fluctuations on diffusive transport in spatially non- 
uniform systems. 
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